/*
 * import_data.cpp
 *
 *  Created on: May 27, 2011
 *      Author: nzhao
 */

#include "special/import_data.h"


double 	* RELATIVE_GAMMA_N;

double	* RELAX_COORD;
double	** HF_MAT_FIRSTPRINCIPLE;
int		N_ATOM_FIRSTPRINCIPLE = 510;

void import_data()
{
	import_coordinates();

    if(IsCore)
    	read_FP_data();
}

void import_coordinates()
{
	int i;
	char comment[200], atom_type[10];
	char full_input_path [300];

	strcpy(full_input_path, INPUT_PATH);
	strcat(full_input_path, INPUT_FILE_NAME);
	strcat(full_input_path, ".xyz");
	ifstream coord(full_input_path);
	if(!coord) assert(0);

	coord >> NBATH;
	COORDINATES		= new double [NBATH * 3];
	ATOM_SPIN		= new int [NBATH];
	ABSOLUTE_GAMMA_N = new double [NBATH];
	RELATIVE_GAMMA_N = new double [NBATH];

	coord >> comment;
	for (i = 0; i < NBATH; i++)
	{
		coord >> atom_type;
		coord >> COORDINATES[ 3 * i + 0];
		coord >> COORDINATES[ 3 * i + 1];
		coord >> COORDINATES[ 3 * i + 2];

		GetProperties(atom_type, i);
		RELATIVE_GAMMA_N[i] = ABSOLUTE_GAMMA_N[i] / GAMMA_N;
	}

	coord.close();
}
void  GetProperties(char * element, int i)
{
    int k, isFound = 0;
    for (k = 0; k < MR_DATA_LENGTH; k++)
    {
        if(strcmp(MR_DATA[k].element, element) == 0)
        {
        	ATOM_SPIN[i] 		= MR_DATA[k].spin_dim;
            ABSOLUTE_GAMMA_N[i] = MR_DATA[k].gamma;
            isFound = 1;
        }

    }

    if(!isFound)
    {
    	cout << element << endl;
    	cout <<"ERROR: missing elements!" << endl;
    	assert(0);
    }
}

void read_FP_data()
{
    int i, j;
    //////////////////////////////////////////////////////////////////////////

    char relaxCoord_file[200];
    strcpy(relaxCoord_file, WORKSHOP);
    strcat(relaxCoord_file, "src/data/RelaxCoord.dat");
    cout << relaxCoord_file << endl;

    RELAX_COORD = new double [3 * N_ATOM_FIRSTPRINCIPLE];

    ifstream relax(relaxCoord_file);
    if(!relax) assert(0);

    for (i = 0; i < N_ATOM_FIRSTPRINCIPLE; i++)
    {
            relax >> RELAX_COORD[ 3 * i + 0 ];
            relax >> RELAX_COORD[ 3 * i + 1 ];
            relax >> RELAX_COORD[ 3 * i + 2 ];
    }

    relax.close();

    //////////////////////////////////////////////////////////////////////////


    char hf_info_file[200];
    strcpy(hf_info_file, WORKSHOP);
    strcat(hf_info_file, "src/FPdata/fhCoeff.dat");
    cout << hf_info_file << endl;

    HF_MAT_FIRSTPRINCIPLE = new double *[N_ATOM_FIRSTPRINCIPLE];

    ifstream hfstream(hf_info_file);
    if(!hfstream) assert(0);

    for (i = 0; i < N_ATOM_FIRSTPRINCIPLE; i++)
    {
            HF_MAT_FIRSTPRINCIPLE[i] = new double [9];
            hfstream >> HF_MAT_FIRSTPRINCIPLE[i][ 0 ];
            hfstream >> HF_MAT_FIRSTPRINCIPLE[i][ 1 ];
            hfstream >> HF_MAT_FIRSTPRINCIPLE[i][ 2 ];

            hfstream >> HF_MAT_FIRSTPRINCIPLE[i][ 3 ];
            hfstream >> HF_MAT_FIRSTPRINCIPLE[i][ 4 ];
            hfstream >> HF_MAT_FIRSTPRINCIPLE[i][ 5 ];

            hfstream >> HF_MAT_FIRSTPRINCIPLE[i][ 6 ];
            hfstream >> HF_MAT_FIRSTPRINCIPLE[i][ 7 ];
            hfstream >> HF_MAT_FIRSTPRINCIPLE[i][ 8 ];
            for(j = 0; j < 9; j++)
            	HF_MAT_FIRSTPRINCIPLE[i][j] *= 0.5;
    }

    hfstream.close();
}

int isNear(int idx)
{
	int i, res;

	double diff, x, y, z;

	x = COORDINATES[3*idx]; y = COORDINATES[3*idx+1]; z = COORDINATES[3*idx+2];

	res = -1;

	for (i = 0; i < N_ATOM_FIRSTPRINCIPLE; i++)
	{
		diff = (RELAX_COORD[i*3] - x)*(RELAX_COORD[i*3] - x) + (RELAX_COORD[i*3+1] - y)*(RELAX_COORD[i*3+1] - y) + (RELAX_COORD[i*3+2] - z)*(RELAX_COORD[i*3+2] - z);
		if (diff < 0.04)
		{
			res = i;
			break;
		}
	}

	return res;
}
